Thermalization Breakdown and Conductivity Improvement within the Interacting 
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Based on the framework of Kubo formulism, we develop the minimally entangled typical ther- 
mal state algorithm to study the temperature and time dependence of current-current correlation 
function in one-dimensional spinless fermion model, taking into account both the electron-electron 
(e-e) intersite interaction and the dynamic disorder induced by classical phonons. Without e-e in- 
teraction, the numerical results, showing an exponential decay of the time dependent correlation, 
could be precisely compared with that from the analytical derivation, namely, from the generalized 
Langevin equation. More importantly, when a strong enough e-e interaction is presence, we find 
a long-time correlation in the regime of small dynamic disorder, indicating the breakdown of ther- 
mal relaxation, which is a typical many-body effect. On the basis of this finding, we show that it 
might be applied to understand the metalliclike charge transport and the abnormal improvement 
of the conductivity with respect to the redoping experiment in K3C60, an organic superconducting 
material. 

PACS numbers: 72.80.Le, 71.10.Fd, 72.10.-d 



Recent progresses of organic superconducting 
materials, such as potassium-doped picene,[l], 2] 
phenanthrene,Q coronene,[iJ and dibenzpentacene, 5] 
has opened a new research subfield in organic electronics, 
due to the following two critical points: The electron- 
intramolecular- vibration interaction, 0, 0] together with 
the intercalant and intermolecular phonons, arc 
substantially respondence to the superconductivity; 
The materials are typical strongly-correlated electron 
systems. [I0I - H2I Both the two statements are from 
first-principle calculations, and an in-depth model 
computation is not found. On the other hand, even 
under room temperature, once doped with alkali metal, 
those originally semiconducting materials become to 
behave metalliclike conductivity, d EI Hi Intuitively, 
the doping of alkali metal has modified the 7r-electron 
structure of the organic molecules, [l[ and the electron- 
electron (e-e) Coulomb interaction becomes to improve 
the conductivity. But in a classical manner, the e-e 
interaction always plays a negative role (blocking) in the 
charge transport. 15| This contradiction implies that, a 
completely quantum description should be addressed for 
this subject, which is the main motivation of this work. 

Beyond (semi-)classical treatment, there have been 
many quantum theories for the transport in organic 
solids. (I6M23I . [25| Most of the works paid attention to 
the subjects, such as the unified bandlike and hopping 
mobility, the static and dynamic disorder, 20 23] 

and the mixed quantum and classical problems. [251] In 
particular, the dynamic disorder model, [20| based upon 
the one-dimensional Holstein-Peierls Hamiltonian with 
both intra- and inter-molecular phonons treated clas- 
sically, was extensively used to comprehensively un- 
derstand the behavior of transport in organic crystals. 



Originally, one used the Ehrenfest method to simulate 
the diffusion behavior of an initially localized electron 
wavepacket [20M22I ] and successfully given the basic car- 
rier's bandlike mobility. [2l| However, as those works were 
mainly working within one-particle picture, they make no 
sense of the fluctuation of the particle number, so that the 
dynamical response could not be evaluated. To fix this 
problem, the Kubo formulism was taken into account, 
in which the mobility is directl y re lated to the current- 
current correlation function. (l8l. |23| It was obtained that, 
both bandlike and hopping transport could be described, 
that is, the localization length decreased quickly when 
the electron-phonon (e-p) coupling increases. [23j At the 
mean time, the e-e interaction is also studied on the mean 
field level. [19( In addition, the ID Holstein model was also 
applied to study the organic superconductivity, since the 
e-p coupling is recognized to be mostly relevant. [26j In 
all, it seems to say that, one can just straightforwardly 
follow this line to study the e-e correlation more com- 
prehensively, which should be the essential character of 
organic superconducting materials. Yet, as we will show 
in this Letter, the breakdown of thermalization 27- 2^| 
induced by the e-e correlation makes the problem quite 
novel. 

The density matrix renormalization group (DMRG), a 
well-known numerical method, is one of the most power- 
ful methods to deal with the one-dimensional strongly- 
correlated systems. 30] In the last decade, lots of ef- 
fort have been put into extending the method to finite- 
temperature problems. [3 u - l35| Advantages from White, 
who is the inventor of DMRG, are made by introduc- 
ing the language of matrix product state and quan- 
tum Monte Carlo method, and a so-called minimally en- 
tangled typical thermal state (METTS) algorithm was 
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established. 35] This new method is highly efficient to 



calculate the thermal quantities of the system of one- 
dimensional spin (and thus spinless fermion) lattices, es- 
pecially under high temperature. Furthermore, if the 
imaginary-time evolution operators in METTS algorithm 
are replaced by its real-time counterpart, which is quite 
straightforward, the method is then applicable for both 
temperature and time dependent problems. This means 
it finally becomes possible to study the thermodynam- 
ics of an organic electronic systems, such as Kubo for- 
mula and time-dependent current-current correlation, 
with both e-p and e-e interaction presence. 

The model we are dealing with is a one-dimensional 
spinless model with near-neighboring e-e interactions, 
and the Holstein e-p coupling, which is treated classically, 
is also taken into account to act as a dynamic disorder. 
The Hamiltonian writes, 



H = H e \ + H ph . 



(1) 



The electronic part is 



H e i = -to^^^Cj +h.c.) + g^up 



TljUj+l, 



(2) 



where ct(c.,) creates (annihilates) an electron on the j-th 

site, Uj the displacement of the j-th site, fij(= ejej) the 
number operator of the electron, to the transfer integral, 
g the e-p coupling constant, and V the intersite e-e inter- 
action. The phonon part of Hamiltonian (fTJ is described 

as 



K 



ph 



(3) 



where K is the elastic constant between neighbor sites, 
and M the mass of a site. All the parameters in the 
model could be determined by first-principle calculations. 
For example, the intermolecular transfer integral to in 
doped corenene is around 30meV, @ and the characteris- 
tic phonon frequency in doped picene is around 18meV. [|| 
But in this work, we will take the dimensionless param- 
eters, that is, t, K and M are all set to unity, such that 
the frequency of phonons is uj = y/K/M = 1. g and 
V will be the main adjustable parameters, and the main 
results are calculated in an open chain with 60 sites. To 
diminish the influence of open boundary condition, we 
have made the to exponential decay on several bonds near 
the boundary. Actually, the ID spinless fermion model 
with e-p coupling has been long-termly studied to under- 
stand the superconductivity in doped fullerenes.[26j It 
was found that, the phase transition from Luttinger liq- 
uid phase to charge density wave occurs around g = 2to 
when t > 1. This critical phenomenon is also found in our 



calculations, but the main physical results in this work 
are within the Luttinger liquid phase. 

Based on the Hamiltonian, we are going to calculate 
the zero-frequency Kubo formula [36l) defined as 



1 



knT 



dt(I(t)I(0)), 



(4) 



where a is the conductivity, T the temperature, I(t) 
(= e iHt/hj e -zHt/h^ foe current operator, and (A) the 
thermal average defined as (A) = Tr(e~ H / klsT A)/ Z with 
Z the partition function. Herein, to define the current 
operator I, one should first define the polarization oper- 
ator P of the center of mass of the electrons, namely, 



(5) 



with e the charge of electron and Rj the position of each 
site. Then I could be defined by 

dP 1 ~ . eatn -e-^ , i , . 

/ =^ = ^] = ^BsV 1 -h.c.), (6) 

3 

with a the lattice constant. 

The current-current correlation function in the Kubo 
formula ([4]) is both time and temperature dependent. 
In order to evaluate it, we apply the METTS algo- 
rithm combined with time-dependent DMRG method 

(tDMRG),^3 

whose basic procedure is as follows. 
Firstly, one initializes a configuration of the occupation 
of each site arbitrarily and then produces a so-called clas- 
sical product state (CPS) as 



|n) = \ni,n 2 , - ■ -rij ■ ••), 



(7) 



with \rij) the local basis on site j. Secondly, the 
imaginary-time evolution operator is acting on the CPS, 
namely, 



l )=P(n)- 1 / 2 e-^\n), 



(8) 



where \4> n ) is the so-called typical thermal state, P{n){= 
(n|e _/3 " f/ \n)) the statistical probability of the state \n) in 
the ensemble, and j3 the inverse of ksT. In a Monte 
Carlo manner, the above steps are iterated, and then we 
get lots of (but still much fewer than the total number of 
state) \4> n ) samplings in hand. Based on these samplings, 
we calculate the expectation of any operator A as 



(9) 



Obviously, when the number of site becomes large, the 
Hilbert space is enlarged drastically. Hence, one should 
follow the standard procedure of tDMRG, that is, the 
evolution operator must be decomposed onto individual 
bonds, and the state must be truncated while scanning 
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as the usual procedure in DMRG method. Meanwhile, 
in White's treatment, to ensure the ergodic hypothesis, 
one should choose another CPS by collapsing \<f) n ) to (the 
arbitrary basis of) each site and iterate the above steps 
for sufficient times. 

Once the thermal equilibrium is reached, that is, 
the typical thermal states \<p n {0)) in the equilibrium 
are obtained, one can easily calculate the time evolu- 
tion by changing the evolution operator with imaginary 
time in to real time. For the aim of the present 
work, we are focusing on computing the time dependent 
current-current correlation function in Kubo formula, 
i.e., {(j> n {Q)\e iHt ' n ie- im l n i\<f> n (<$)). Hence, one should 
first make the current operator / acting onto \<fi n (0)), and 
then calculate the time evolution of the obtained state. 
At each time step t, one again act I on the new state 
and calculate its overlap with the state e~ lHt / h \4> n (Qi)) . 
Since / could be decomposed onto each bond of the lat- 
tice, naively, one can just act it on the individual bond in 
each scanning step and sum up the targeting states. But 
this idea is quite inefficient for tDMRG, because we need 
to target too many states. Our treatment is on the basis 
of assumption of translating invariance, namely, at the 
initial moment, we only act the current operator on the 
central bond of |0„(O)), which should be the most pre- 
cise bond in the approximation of tDMRG. Based on this 
simplicity, only two states \4> n {0)) and I c |0n(O)) (with I c 
the current operator on the central bond) are necessary 
to be targeted, which makes the procedure much more 
efficient. 

Up to now, the remaining thing is to treat the phonon 
part of the system. We follow the usual procedure of 
Ehrenfest method, which is the standard method for cal- 
culating the time evolution of mixed classical-quantum 
systems. That is, we first choose an initial configuration 
of {uj} and {iij} from the Gaussian distribution with 
variance ksT / K and ksT/M, say the equilibrium dis- 
tribution of vibrations. This configuration could be sub- 
stituted into the Hamiltonian of electron part, and one 
follows the common procedure to calculate the Hellman- 
Feynman force that electrons act on the sites, such that, 
the influence of phonon comes into the theory. 

Here, we need to say that, since the current opera- 
tor / docs not change the total number of electron and 
the spin degree of freedom is ignored, the precision of 
the present numerical calculation should be at least in 
the same order with the previous tDMRG works of us in 
the Hubbard chain and the conjugated polymer. [38[ How- 
ever, the evolution is now for the typical thermal states 
rather than the ground state. Together with the numer- 
ical errors from the Ehrenfest method for the phonon 
part, the accumulation of error is very fast. Even if we 
decrease the number of truncated states in METTS, it 
seems no significant improvement is obtained, so that 
more efforts should be devoted to the method itself. But 
still, within short time scale, many interesting results 
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FIG. 1: (a) The current-current correlation function (in the 
unit of (eato) 2 /2N c h 2 ) versus time (in the unit of h/to) with 
various e-p coupling, (b) The averaged velocity from the 
generalized Langevin equation. (3^] The relaxation time for 
g = 0.25 is around 50. 

have been found. 

In Fig. Ufa), we nrs ^ show the results of current- 
current correlation function without e-e interaction under 
ksT = 2. For g — 0, i.e., the disorder is absent, the elec- 
tron is completely free. One can easily prove that, in this 
case, I commutes with H, so the correlation function is 
unchanged with time. As we see that, the curve remains 
a constant up to t = 20, but due to the accumulation of 
numerical error, the curve goes down quickly after that 
point. Hence, t = 20 should be a point of justification 
that, before it the numerical results are credible. In addi- 
tion, we find that, following g increasing, the decay of the 
correlation becomes much faster, and above g — 2 (not 
shown), the correlation acts nearly a sudden quenching, 
which means it is within an insulating phase. 

To check the correctness of the present numerical algo- 
rithm, a method beyond the framework of Monte Carlo 
sampling is needed. Here, we adopt the generalized 
Langevin equation (GLE) to calculate the time depen- 
dence of thermal averaged velocity. Although not ex- 
actly the same, the oscillation and decay of the averaged 
velocity are expected to follow the similar trajectory with 
that of current-current correlation function. To compare 
with the METTS results, in Fig. [TJb), we show the result 
of V(t) from GLE with the same g's, respectively. Qual- 
itatively, both the wavy and decay shape agree with the 
result from METTS. And more importantly, when g de- 
creases by the same times, the decay velocity shows very 
close value, which could not be adjusted by any param- 
eters. These results not only allow us to doubly check 
the correctness of the numerical algorithm, but also pro- 
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FIG. 2: Time dependence of current-current correlation func- 
tion for various g and V. Inset of (c) shows the little influence 
of the finite length scale. 



vide more information that METTS can not give. For 
example, the GLE result shows the relaxation behavior 
for all g after a long time evolution with the relaxation 
time is around 50. This complement is very important 
to understand the following results. 

Now, we are on the stage to show the influence of 
e-e interaction. In Fig. [2l we show the time depen- 
dent correlation for various V and g. Compared with 
V = case, when g is larger than 0.5, the velocity of 
decay increases with increasing V. This is easy to under- 
stand, since the system tend to the insulating phase. But 
quite interestingly, as one can see from (a) to (d) that, 
when g — 0.25, although the velocity of decay still tends 
to increase in the very beginning of evolution, it seems 
that the curve does not vanish but relax to a finite value 
(around 0.15) after a long time. Of course, one may ar- 
gue that, our numerical method can not show the result 
of long-time evolution, but we can still estimate from 
the present result that, the relaxation time for V = 2 
should be at least much longer than that in uncorrelated 
system (s> 50). Meanwhile, to exclude the influence of 
the chain length, in the inset of Fig. [Jfc) , we show the 
result for different site number, and obviously, they are 
almost the same. Hence, this long-time memory effect 
of correlation function is the main finding of this work. 
It is non- asymptotic from the uncorrelated system, since 
it implies the thermalization breakdown induced by the 
e-e interaction. [27]] Actually, the breakdown of thermal 
relaxation in fermion systems was widely studied very 
recently. [H, [29| That is, when the main-body correla- 
tion enters into the low-dimensional fermion system, the 
relaxation to thermal equilibrium should be strongly re- 
lated to the initial state, due to the complex entangle- 
ment between equilibrium and nonequilibrium correla- 
tions in correlated systems. [2^] The present results agree 
with these works and add more information on this sub- 
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FIG. 3: Temperature dependence of conductivity (in the unit 
of to(ea) 2 /2h) calculated from Kubo formula for various V. 



ject. 

Intuitively, since the conductivity is related to the in- 
tegral of current-current correlation function within the 
framework of Kubo formula, the long-time memory effect 
found above should contribute a lot to the conductivity 
and improve it. This statement obviously contradicts to 
the conventional point of view, as the Coulomb repulsion 
in the material with small disorder is always recognized 
to act as a blocking and thus will lead to the reduction of 
conductivity. [l5| Whereas, the improvement of conduc- 
tivity was indeed found in K3C60 by redoping of alkali 
metals, as shown in (ljj . The redoping process is, in our 
opinion, to increase the concentration of alkali metallic 
atoms and thus the Coulomb interaction among them. 
The improvement of conductivity is surprising and full 
of meaning, since it is closely related to the increase of 
Tc of these organic superconducting materials. 13j On 
the basis of our present finding of the memory effect, 
we could then provide a possible explanation on the ex- 
periment. In Fig. [31 we show the main results of the 
temperature dependence of conductivity. Here, due to 
the limit of numerical method, the time integral in the 
Kubo formula ((4]) is only within t < 20. This is equiva- 
lent to considering an additional static disorder l/r s , i.e., 
the time integral in Kubo formula becomes [3| 



dt -> / dte~ (t/T 



(10) 



It is found that, when V < 2, the result is quite the same 
with the bandlike behavior in common dynamic disor- 
der model. [2(J When V > 2, the relationship between 
conductivity and temperature becomes to metalliclike, 
namely, a ~ 1/T. More importantly, the curve of con- 
ductivity stops decreasing but becomes close with each 
other. Even in the high temperature regime and V = 4, 
we find a cross of the curves and improvement of the con- 
ductivity, which are consistent with that of doped and 
annealed cases as shown in the Fig. 3 of [l3j]. Therefore, 
we could now conclude that, the thermalization break- 
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down induced by the e-e interaction matters in the charge 
transport of alkali-metal-doped organic materials. 

Finally, we would like to briefly discuss the supercon- 
ductivity in organic solids. Different from their inor- 
ganic counterparts, organic materials have much more 
diverse vibrational modes, such that the decoherence pro- 
cess could easily kill the phase correlation of the electron 
wavefunctions and thus the tendency of superconductiv- 
ity. To avoid this effect, the high frequency part of the 
phonons, which is the main source of decoherence, 24, 4pj 
must be largely suppressed. The present work provides a 
possibility, that is, the long-time memory of electric cur- 
rent induced by e-e interaction should be a shield against 
the decoherence and contribution to the conduction of 
these materials. 

In summary, we have used the METTS algorithm to 
calculate the temperature and time dependent current- 
current correlation function in a one-dimensional spin- 
less model with both e-e and e-p interaction taken into 
account. Via the comparison with analytic results, we 
state that, this very new method works well within a 
short time scale. Then we study the influence of e-e in- 
teraction and find a long-time memory effect, say, the 
thermalization breakdown of the system. Based on this 
finding, we calculate the temperature dependent mobility 
and show that, following the increase of e-e interaction, 
the mobility will behave a slight enhancement, which was 
also found in the experiment. 
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In order to make sense of the long-time behavior for V = 0, we adopt the generalized Langevin equation (GLE)ji^ 
The basic idea is to separate the electronic Hamiltonian into the center-of-mass and relative part, namely, 

pi 

6 k q 



where P is the momentum of the center of mass, R the conjugate variable of P, M e = Nm e with N the number of 
electron and m e the mass of a single electron, the energy with k the wave vector of "relative" electrons, ct(cfc) the 
creation (annihilation) operator of relative electrons with p q (= ^2k^k+ q ^k) the corresponding density operator, and 
q the wave number. It could be demonstrated that, the operators of center-of-mass and relative part commute with 
each other^ Meanwhile, only the contribution of small q is assumed to be dominant, so that the relative part acts as 
a source of fluctuations for the center of mass via the assistance of phonons. 

Within this Hamiltonian, we are going to evaluate the derivative of V{= P/M e ) at time t. In the interaction 
picture, it could be written as (H = l)i 

^ = F{t)- j\t'[F{t),H ep (t% (2) 

where H ep is the last term of ([1}, an d F(t) is the fluctuating force operator defined as 

F(i) = -i^e"'V«,f,, (3) 
q' 

Since q is small, we can expand the exponential term to the linear order, that is, exp(zgi?) = 1 + iqR. Then Eq. ([5]) 
becomes 

dV rt 
~dt 



= F(t)- dt'g\R{t) + R{t'))Y,Q 2 u%t').[p q {t)rpS% (4) 
Jo q 

where the zero-order term is absorbed into the fluctuating force, and the high-order term of q is neglected, q' is equal 
to q in order to ensure the conservation of the total momentum. 

Following the Ehrenfest theorem, the dynamics of classical phonons could be treated as u q (t,t') ~ ksT sin(u;(i — 
t'))/K with lo its frequency and K the elastic constant^ The last term, say [p q (t), p q (t')], behaves as a relaxation 
term and is significant when t = t'. It is not easy to derive an explicit expression for this term, but it could be 
estimated to be ~ n/q 2 with n the density of electron^ In all, we then define a dimensionless parameter £, which 
denotes the thermal average of the summation, i.e.,(^ q 2 u q (t, t') ■ [p q (t), p q (t')]) = £ sin(w(< — t')). Actually, £ equals 

to the intensity of the fluctuating force and is proportional to hsT^^ Finally, by taking the thermal average and 
using integration by parts, Eq. ^ becomes the GLE as 

^ = Ht) - £g 2 [2R(t) - J* dt' cos(uj(t - t'))V(t% (5) 

where we have used Vif) = dR(t) / dt, and the cosine term is absorbed into the fluctuating force. Here £ is an adjustable 
parameter, and when it is sufficiently small, R(t) ~ V(t)t. To compare with the numerical results, we will set £ to be 
0.04. Then we take the thermal average for all the quantities and make (F) vanishing. The numerical integration is 
now available for V(t) when we consider V^(0) = 1 as the initial condition. 
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